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Abstract 

We present a simple Markov model of spiking neural dynamics that can be analytically solved to 
characterize the stochastic dynamics of a finite-size spiking neural network. We give closed-form 
estimates for the equilibrium distribution, mean rate, variance and autocorrelation function of the 
network activity. The model is applicable to any network where the probability of firing of a neuron 
in the network only depends on the number of neurons that fired in a previous temporal epoch. 
Networks with statistically homogeneous connectivity and membrane and synaptic time constants 
that are not excessively long could satisfy these conditions. Our model completely accounts for 
the size of the network and correlations in the firing activity. It also allows us to examine how the 
network dynamics can deviate from mean-field theory. We show that the model and solutions are 
applicable to spiking neural networks in biophysically plausible parameter regimes. 



* Neural Computation, in press 



1 



INTRODUCTION 



Neurons in the cortex, while exhibiting signs of synchrony in certain states and 

tasks QQ0.Q 

, mostly fire stochastically or asynchronously [5J . Previous theoretical and 
computational work on the stochastic dynamics of neuronal networks have mostly focused 
on the behavior of networks in the infinite size limit with and without the presence of exter- 
nal noise. The bulk of these studies utilize a mean-field theory approach, which presumes 
self-consistency between the input to a given neuron from the collective dynamics of the 
network with the output of that neuron and either ignores fluctuations or assumes that the 
ductearions obey a prescribed statists distribution, e.g. fl fl 0, fl U Q aad Q for 
a review. Within the mean-field framework, the statistical distribution of fluctuations may 
be directly computed using a Fokker-Planck approach [13I 
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estimated from the response of one individual neuron submitted to noise 

Although, correlations in mean-field or near mean-field networks can be nontrivial 
S, H], in general, mean-field theory is strictly valid only if the size of the network is large 
enough, the connections are sparse enough or the external noise is large enough to decorrelate 
neurons or suppress fluctuations. Hence, mean-field theory may not capture correlations in 
the firing activity of the network that could be important for the transmission and coding 
of information. 



It is well known that finite-size effects can contribute to fluctuations and correlations [28 ]. 
It could be possible that for some areas of the brain, small regions are statistically homo- 
geneous in that the probability of firing of a given neuron is mostly affected by a common 
input and the neural activity within a given neighborhood. These neighborhoods may then 
be influenced by finite-size effects. However, the effect of finite size on neural circuit dy- 
namics is not fully understood. Finite-size effects have been considered previously using 



expansions around mean-field theory 
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29j. It would be useful to develop analytical 



methods that could account for correlated fluctuations due to the finite size of the network 
away from the mean- field limit. 

In general, finite-size effects are difficult to analyze. We circumvent some of the diffi- 
culties by considering a simple Markov model where the firing activity of a neuron at a 
given temporal epoch only depends on the activity of the network at a previous epoch. 
This simplification, which presumes statistical homogeneity in the network, allows us to 
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produce analytical expressions for the equilibrium distribution, mean, variance, and auto- 
correlation function of the network activity (firing rate) for arbitrary network size. We find 
that mean-field theory can be used to estimate the mean activity but not the variance and 
autocorrelation function. Our model can describe the stochastic dynamics of spiking neural 
networks in biophysically reasonable parameter regimes. We apply our formalism to three 
different spiking neuronal models. 

THE MODEL 

We consider a simple Markov model of spiking neurons. We assume that the network 
activity, which we define as the number of neurons firing at a given time, is characterized 
entirely by the network activity in a previous temporal epoch. In particular, the number 
of neurons that fire between t and t + At only depends on the number that fired between 
t — At and t. This amounts to assuming that the neurons are statistically homogeneous 
in that the probability that any neuron will fire is uniform for all neurons in the network 
and have a short memory of earlier network activity. Statistical homogeneity could arise for 
example, if the network receives common input and the recurrent connectivity is effectively 
homogeneous. We note that our model does not necessarily require that the network ar- 
chitecture be homogeneous, only that the firing statistics of each neuron be homogeneous. 
As we will show later, the model can, in some circumstances, tolerate random connection 
heterogeneity. Short neuron memory can arise if the membrane and synaptic time constants 
are not excessively long. 

The crucial element for our formalism is the gain or response function of a neuron p(i, t), 
which gives the probability of firing of one neuron within t and t + At given that % neurons 
have fired in the previous epoch. The time dependence can reflect the time dependence of 
an external input or a network process such as adaptation or synaptic depression. Without 
loss of generality, we can rescale time so that At = \. 

Assuming p(i,t) is the same for all neurons, then the probability that exactly j neurons 
in the network will fire between t and t + 1 given that i neurons fired between t — 1 and t is 

P (x(t + 1) = j\x(t) =i) = c j NP {i, ty(i - P (i, t)) N -*, (i) 

where X(t) is the number of neurons firing between t — 1 and t and C% is the binomial 
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coefficient. Equation ([I]) describes a Markov process on the finite set {0, ..,N} where N is 
the maximum number of neurons allowed to fire on a time interval At. The process can be 
re-expressed in terms of a time dependent probability density function (PDF) of the network 
activity P(t) and a Markov transition matrix (MTM) defined by 

M^t) = p (X(t + l) =j\X(t) = i) = C k N p(i, ty(l-p(i,t)) N -i. (2) 

P(t) is a row vector on the space [0, 1]^ that obeys 

P(t + 1) = P(t)M(t) (3) 

For time invariant transition probabilities, p(i,t) = p(i) and for fixed N, we can write 

P(t) = P(0)M\ (4) 

Our formalism is a simplified variation of previous statistical approaches, which use renewal 



models for neuronal firing 



27 



3l|. In those approaches, the probability for a given 
neuron to fire depends on the refractory dynamics and inputs received over the time interval 
since the last time that particular neuron fired. Our model assumes that the neurons are 
completely indistinguishable so the only relevant measure of the network dynamics is the 
number of neurons active within a temporal epoch. Thus the probability of a neuron to 
fire depends only on the network activity in a previous temporal epoch. This simplifying 
assumption allows us to compute the PDF of the network activity, the first two moments, 
and the autocorrelation function analytically. 



TIME INDEPENDENT MARKOV MODEL 



We consider a network of size TV with a time independent response function, p(i, t) = p(i), 
so that (jl]) specifies the temporal evolution of the PDF of the network activity. If < p(i) < 1 
for all i G [0, N] (i.e. the probability of firing is never zero nor one), then the MTM is a 
positive matrix (i.e. it has only non-zero entries). The sum over a row of the MTM is 
unity by construction. Thus the maximum row sum matrix norm of the MTM is unity 
implying that the spectral radius is also unity. Hence, by the Frobenius-Perron Theorem, 
the maximum eigenvalue is unity and it is unique. This implies that for all initial states of 
the PDF, P(t) converges to a unique left eigenvector [i called the invariant measure, which 
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satisfies the relation 

H = fiM. (5) 

The invariant measure is the attractor of the network dynamics and the equilibrium PDF. 
Given that the spectral radius of M is unity, convergence to the equilibrium will be expo- 
nential at a rate governed by the penultimate eigenvalue. We note that if the MTM is not 
positive then there may not be an invariant measure. For example if the probability of just 
one possible transition is zero, then the PDF may never settle to an equilibrium. 

The PDF specifies all the information of the system and can be found by solving (j3j). 
From , we see that the invariant measure is given by the column average over any arbitrary 
PDF of lim^oo M t . Hence, the infinite power of the MTM must be a matrix with equal rows, 
each row being the invariant measure. Thus, a way to approximate the invariant measure to 
arbitrary precision is to take one row of a large power of the MTM. The higher the power, 
the more precise the approximation. Since the convergence to equilibrium is exponential, a 
very accurate approximation can be obtained easily. 

Mean and variance 

In general, of most interest are the first two moments of the PDF so we derive expressions 
for the expectation value and variance of any function of the network activity. Let X 6 
[0, ...N] be a random variable representing the network activity and / be a real valued 
function of X. The expectation value and variance of / at time t is thus 

N 

(/(*)>* = £/(*W*) (6) 

k=0 

and 

Var t (f(X)) = (f(Xr) t -(f(X))l (7) 

where we denote the kth element of vector P(t) by Pk{t). We note that in numerical 
simulations, we will implicitly assume ergodicity at equilibrium. That is, we will consider 
that the expectation value over all the possible outcomes to be equal to the expectation over 
time. 

Inserting ([2j) into (jHj) and using the definition of the mean of a binomial distribution, the 
mean activity (X) t has the form 

(X) t = N(p(X)) t ^. (8) 
5 



The mean firing rate of the network is (X) t /At. We can also show that the variance is given 
by 

Var t (X) = N(p(l - p)) t -x + N^Var^ip). (9) 

Details of the calculations for the mean and variance are in Appendix . Thus the mean 
and variance of the network activity are expressible in terms of the mean and variance of 
the response function. At equilibrium, ([H]) and Q are (A) M = N(p) fM and Var^X) = 
N(p(l — p))^ + ^Var^p), respectively. Given an MTM, we can calculate the invariant 
measure and then from that obtain all the moments. 

If we suppose the response function to be linear, we can compute the mean and variance 
in closed-form. Consider the linear response function 

p(*)=Po + ^f^ (10) 

for X G [0, N]. Here po G [0, 1] is the probability of firing for no inputs and (q — Po)/Nq is 
the slope of the response function. Inserting (ITUl) into (JED, gives 

W, = iVp„+^(I), (11) 

Solving for (X)^ gives the mean activity 

(X), = Nq, (12) 

Substituting (TTDT) into ([9]) leads to the variance 

Va '^ = N l -tit IN < 13 > 

where A = (q — Po)/q- The details of these calculations are given in Appendix . 

The expressions for the mean and variance give several insights into the model. From 
(flOl) . we see that the mean activity (fT2l) satisfies the condition 



p(Nq) = q. (14) 

Hence, on a graph of the response function versus the input, the mean response probability 
is given by the intersection of the response function and a line of slope 1/N, which we denote 
the diagonal. Using ffT2l . f|T4|) can be re-expressed as 

(X), = Np((X)^ (15) 
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In equilibrium, the mean response of a neuron to (X) neurons firing is (X)/N. Hence, 
for a linear response function, the mean network activity obeys the self-consistency condi- 
tion of mean-field theory. This can be expected because of the assumption of statistical 
homogeneity. 

The variance fll3p does not match a mean-field theory that assumes uncorrelated statis- 
tically independent firing. In equilibrium, the mean probability to fire of a single neuron in 
the network is q. Thus, each neuron obeys a binomial distribution with variance q(l — q). 
Mean-field theory would predict that the network variance would then be given by Nq(l — q). 
Equation ( Fl3l) shows that the variance exceeds the statistically independent result by a fac- 
tor that is size dependent, except for q = p , which corresponds to an unconnected network. 
Hence, the variance of the network activity cannot be discerned from the firing charac- 
teristics of a single neuron. The network could possess very large correlations while each 
constituent neuron would exhibit uncorrelated firing. 

When N » A 2 , the variance scales as N but for small N, there is a size-dependent 
correction. The coefficient of variation approaches zero as 1/ \^N as N approaches infinity. 
When A = 1, the slope of the response function is 1/N and coincides with the diagonal. 
This is a critical point where the equilibrium state becomes a line attractor 32|. In the limit 
of N — > oo, the variance diverges at the critical point. At criticality, the variance of our 
model has the form N 2 q(l — q), which diverges as the square of N. Away from the critical 
point, in the limit as N goes to infinity, the variance has the form (Nq(l — q))/{l — A 2 ). 
Thus, the deviation from mean-field theory is present for all network sizes and becomes more 
pronounced near criticality. 

The mean-field solution of fTH|) is not a strict fixed point of the dynamics (jlj) per se. For 
example, if the response function is near criticality, possesses discontinuities, or crosses the 
diagonal multiple times then the solution of (fl4|) may not give the correct mean activity. 
Additionally, if the slope of the crossing has a magnitude greater than 1/N, then the crossing 
point would be "locally unstable". Consider, small perturbations of the mean around the 
fixed point given by ( |T2l) 



(X) t = Nq + v(t) (16) 
The response function fflOl) then takes the form 

p(X) = q+±v (17) 
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where A = (q — po)/q. Substituting (|T6|) and (|T7|) into ([8]) gives 



(v)t = Hv) 



t-i 



(18) 



Hence, v will approach zero (i.e. fixed point is stable) only if |A| < 1. 

Finally, we note that in the case that there is only one solution to (TH1) (i.e. response 
function crosses the diagonal only once), if the function is continuous and monotonic then 
it can only cross with A < 1 (slope less than 1/N) since < p < 1. Thus, a continuous 
monotonic increasing response function that crosses the diagonal only once, has a stable 
mean activity given by this crossing. This mean activity corresponds to the mode of the 
invariant measure. 

We can estimate the mean and variance of the activity for a smooth nonlinear response 
function that crosses the diagonal once by linearizing the response function around the 
crossing point (q = p(Nq)). Hence, near the intersection 



Additionally, the crossing point is stable if |A| < 1. These linearized estimates are likely to 
break down near criticality (i.e. A near one). 

We show an example of a response function and resulting MTM in Figured) Starting with 
a trivial density function (all neurons are silent), we show the time evolution of the mean 
and the variance in figure [DC. We show the invariant measure in figure [Tp. The equilibrium 
state is given by the intersection of the response function with the diagonal. We see that 
the mode of the invariant measure is aligned with the mean activity given by the crossing 
point. 

Autocovariance and autocorrelation functions 




(19) 



where A = Njj^(Nq). Using flTJ) and ([13]) then gives 



(X), = Nq 
Var,{X) = -1 



1 - A 2 + A 2 /A 



Nq(l - q) 



(20) 
(21) 



We can also compute the autocovariance function of the network activity: 
Cov(r) = <(X(t) - (X) M )(X(t + r) - (X),)) = (X(t)X(t + r)> - (X 2 ) 
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(22) 



Noting that 



(X(t)X(t + r)) = X; Y,3 k P(X(t + 7") = = J>(*(*) = 3) 



(23) 



i ft 



where p(X(t + r) = fc|X(t) = j) = ML., we show in Appendix , that at equilibrium 



where A is the slope factor (N times the slope) of the response function evaluated at the 
crossing point with the diagonal. The autocorrelation function AC(t) is simply A T . 

The correlation time is given by 1/ In A. At the critical point (A = 1), the correlation time 
becomes infinite. For a totally decorrelated network, A = 0, giving a Kronecker delta function 
for the autocorrelation as expected. For an inhibitory network (A < 0), the autocorrelation 
exhibits decaying oscillations with period equal to the time step of the Markov approximation 
(i.e. Cov(l) < 0). (Since we have assumed p > 0, an inhibitory network in this formalism is 
presumed to be driven by an external current and the probability of firing decreases when 
other neurons in the local network fire.) We stress that these correlations are only apparent 
at the network level. Unless the rate is very high, a single neuron will have a very short 
correlation time. 

Multiple Crossing Points 

For dynamics where < p < 1, if the response function is continuous and crosses the 
diagonal more than once, then the number of crossings must be odd. Consider the case 
with three crossings with p(0) < p(N). If the slopes of all the crossings are positive (i.e. a 
sigmoidal response function), then the slope factor A of the middle crossing will be greater 
than one while the other two crossings will have A less than one implying two stable crossing 
points and one unstable one in the middle. Figure [2A and B provide an example of such 
a response function and its associated MTM. The invariant measure is shown in figure [2C. 
We see that it is bimodal with local peaks at the two crossing points with A less than one. 
The invariant measure is well approximated by any row of a large power of the MTM (figure 
I2D). If the stable crossing points are separated enough (beyond a few standard deviations), 
then we can estimate the network mean and variance by combining the two local mean and 
variances. Assuming that the stable crossing points are q± and q 2 (q% < q 2 ) with slope factors 



Cov(r) = A T Var M (X) 



(24) 
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Ai and A2. Then the mean and the variance are: 

= (25) 

Var ^ X) = 2 [l-M + M/N + 1 - AI + XV N ) + N [—2-) (26) 

FAST LEAK SPIKING MODEL 

We now consider a simple neuron model where the firing state obeys 

1 — 

Vi(t) = H(I(t) + JikV k (t - 1) + S(t) - 6) (27) 
iN k=i 

where H(-) is the Heaviside step function, I(t) is an input current, S(t) = N(0, a) is uncor- 
related Gaussian noise, 9 is a threshold, and is the connection strength between neurons 
in the network. When the combined input of the neuron is below threshold, Vi — (neuron 



is silent) and when the input exceeds thresho 
version of a spike response or renewal model 



d, Vj = 1 (neuron is active) . This is a simple 
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si 



In order to apply our Markov formalism, we need to construct the response function. 
We assume first that the connection matrix has constant entries J. We will also consider 
disordered connectivity later. The response function is the probability that a neuron will 
fire given that n neurons fired previously. The response function at t for the neuron is then 
given by 



00 

p(n,t) = -f= L r , , , /w e~ x /2 dx (28) 



B-I(t)-nJ/N 



If we assume a time independent input then we can construct the MTM using ()2]). The 
equilibrium PDF is given by (j3J), and can be approximated by one row of a large power of 
the MTM. 

Mean, variance and autocorrelation 

Imposing the self-consistency condition ( fl4l) on the mean activity Nq, where q = 
and p is given by (128"]) leads to 

? = 4= lZ- gJ e- x2/2d *- ( 29 ) 



'2ti 
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Taking the derivative of p(n, t) in (|28|) at n = Nq gives 



J (B-I-Jq) 2 

A = — =e '. (30) 

cr\/27r 



Using this in (|2ip gives the variance. 

Equation (1301) shows the influence of noise and recurrent excitation on the slope factor 
and hence variance of the network activity. For A < 1, the variance fl2T]) increases with A. 
An increase in excitation J increases A for large J but decreases A for small J. Similarly, a 
decrease in noise a increases A for large a but decreases it for small a. Thus, variability in 
the system can increase or decrease with synaptic excitation and external noise depending 
on the regime. The autocovariance function, from (124")) is 

Bifurcation diagram 

Bifurcations occur when the stability of the mean-field solutions (i.e. crossing points) for 
the mean activity changes. We can construct the bifurcation diagram on the I — J plane by 
locating crossing points where A = 1. We first consider rescaled parameters X = (6 — I) ja 
and J = J I a. Then f l2"$j) becomes 

q = ^ / e~ x /2 dx (32) 



and (13"0"j) becomes 



X= J-^ (33 ) 

V Z7T 



Solving (!33l) for g and inserting into (|32|) gives 



* = -i Tr-Tr-^ 2 / 2 ^ ± \ H^-) (34) 



'toJ± y / H &) V v 2vrA 

The bifurcation points are obtained by setting A = 1 in fl34|) and evaluating the integral 



The solutions have two branches that satisfy the conditions J > v 2n and X > y 7r/2. Figure 
[3JA. shows the two dimensional bifurcation diagram on the X — J plane. The intersection of 
the two branches at J = y/2n and X = */| is a co- dimension two cusp bifurcation point (i.e. 



satisfies Np'(x) = 1, Np"(x) = with genericity conditions 33j). Between the two branches 



there are three crossing points (two of which are stable), outside there is one stable crossing 
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point. Each traversal of a branch yields a saddle node bifurcation as seen in figures [3|3 and 
C for the vertical and horizontal traversals (shown in the inset to figure [3]A.) respectively. 
Traversing through the critical point along the diagonal line in the inset of figure [3J\ shows 
a pitchfork bifurcation. We have assumed that J > (i.e excitatory network). By taking 
J — > —J we obtain the same bifurcation diagram for the inhibitory case and the above 
conditions still hold using the absolute value of J. 

We compare our analytical estimates from the Markov model with numerical simulations 
of the fast leak model (12?|) in figure HI Figure H]A. shows that the equilibrium PDF de- 
rived from the MTM matches the PDF generated from the simulation. Figure HB shows 
a comparison of A T versus the simulation autocorrelation function. Figure HP shows the 
mean and variance of the activity versus the connection weight J. The theoretical results 
match the simulation results very well. On the same plot, we show the mean-field theory 
estimate where the neurons are assumed to be completely uncorrelated. For low recurrent 
excitation and hence low activity, the correlations can be neglected and the system behaves 
like an uncorrelated neural network. However, as the excitation and firing activity increase, 
correlations become more significant and the variance grows quicker than the uncorrelated 
result. For very strong connections, the firing activity saturates and the variance drops. 
In that case, mean-field theory is again an adequate approximation. Figure HP shows the 
variance versus N for a network at the critical point, J = \[2tx and X = where A = 1 
at the crossing point. We see that the variance grows approximately as N 3 ^ 2 . This is faster 
than the mean-field result but slower than the prediction for a linear response function. The 
breakdown of the linear estimate is expected near criticality. 

Figure [5] shows the simulation autocorrelation for both the network and one neuron 
(chosen randomly) for iV = 1000 and iV = 10, 000. The parameters were chosen so the 
network was at a critical point (A = 1). We see that the correlation time of an individual 
neuron is near zero and much shorter than that of the network. With increasing N, we find 
that the correlation time increases for the network and decreases for the single neuron. Thus, 
a single neuron can be completely uncorrelated while the network has very long correlation 
times. 
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Bistable regime 



In the regime where there are three crossing points, the invariant measure is bimodal. 



Hence, the network activity is effectived 
posed as a model of working memory 



y 



bistable. Bistable network activity has been pro- 
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371 ] . The two local equilibria are locally 



stable but due to the finite-size fluctuations, there can be spontaneous transitions from one 
state to the other. Figure 0\ shows the activity over time for a bistable network with 
iV = 50. The associated PDF is shown in figure 03. No external input has been added 
to the network. For N = 1000, the switching rate is sufficiently low to be neglected. The 
resulting probability density function for this network is displayed on figure [HP- 

We can estimate the scaling behavior on iV of the state switching rate and hence lifetime 
of a memory by supposing the finite-size induced fluctuations can be approximated by un- 
correlated Gaussian noise. We can then map the switching dynamics to barrier hopping of 
a stochastically forced particle in a double potential well V(x) with a stationary PDF given 
by the invariant measure \i. The switching rate is given by the formula [38J 

_! / (V(a)-V(c)) 



r = (2ir)-yV"(c)\V"(a)\e ~ 2 (35) 

where V = —a 2 ln(/x), \x is the invariant measure, c is a stable crossing point, a is the 
unstable crossing point, and a 2 is the forcing noise variance around c, which we take to be 
proportional to the variance of the network activity. For the general case, we cannot derive 
this rate analytically but we can obtain an estimate by assuming that the invariant measure 
is a sum of two Gaussian PDFs whose means are exactly the two stable crossing points. 



where a 2 = [zf^x^jv f° r ^ = 1, 2. Using this in fl35l) then gives 

r oc {le~ kN , (37) 

for a constant k. The switching rate is an exponentially decreasing function of N. We 
computed the switching rate for iV G [0; 150] for the parameters of the bistable network 
above. The results are displayed on figure 03. For the parameters in 03, a fit shows that k 
is of order 0.05. For these parameters, switching between states will not be important when 
N » 1/k ~ 20. 
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Disordered connectivity 

We now consider the dynamics for disordered (random) connection weights J^- drawn 
from a Gaussian distribution (N(J,a 2 j)). Soula et al. (2006) showed that the input to a 
given neuron for a disordered network can be approximated by stochastic input. We thus 
choose random input drawn from the distribution, N(kJ,kaj), for k neurons having fired 
in the previous epoch. The response function (for time independent input) then obeys 

pin) = / _, e~ x /2 dx. (38) 

Applying the self-consistency condition (fT4"j) for the mean firing probability q gives 

q = -^ - e~ x /2 dx (39) 

The variance is given by (1211) with 

< 1 f . J e - I - ql V-S^ (40) 

Figure [7] shows the activity mean and variance as a function of the connection weight vari- 
ance. There is a close match between the prediction and the results generated from a direct 
numerical simulation with 100 neurons. We note that a possible consequence of disorder in 
a neural network is a spin glass where many competing activity states could co-exist and 
so the resulting activity is strongly dependent on the initial condition [39]. However, we 
believe that the stochastic forcing inherent in our network is large enough to overcome any 
spin glass effects. In the language of statistical mechanics, our system is at a high enough 
temperature to smooth the free energy function. 

INTEGRATE- AND-FIRE MODEL 

We now apply our formalism to a more biophysical integrate-and-fire model. We consider 
two versions. The first treats synaptic input as a current source and the second considers 
synaptic input in the form of a conductance change. We apply our Markov model with a 
discrete time step chosen to be the larger of the membrane or synaptic time constants. 
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Current-based synapse 



We first consider an integrate-and-fire model where the synaptic input is applied as a 
current. The membrane potential obeys 

Devs , „ J N 



T- 



I-V + -J2a(t-t f ) + Z(t) (41) 



dt N , 

V <t 



where a(t) = e~~ jr s if t > and zero otherwise, J is a constant input, Z is an uncorrelated 
zero-mean white noise with variance a, J is the synaptic strength coefficient, N is the 
network size, and v are the firing times of all neurons in the network. These firing times are 
computed whenever the membrane potential crosses a threshold Vg, whereupon the potential 
is reset to Vr s . After firing, neurons are prevented from firing during an absolute refractory 
period r. The parameter values are r = 1 ms, r s = 1 ms,V^. = —65 mV, Vg = —60 mV, 
Vr s = —80 mV and r = 1 ms. 

The network dynamics of this model had been studied previously in detail with a Fokker- 



Planck formalism 
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21| . These analyses showed that if the input current is an uncorrelated 
white noise N(/j,i, a]) and r s < 1 then the response function of the neuron in equilibrium is 
given by 

Vg-llj-JX/N 

r + ^ f v * JXIN due u \l + erf(u)) (42) 



p(X) 

The mean activity is again given by Nq where (p(X))^ = q, and q is a solution of ([T 



« kJl * due" 2 (l + erf(u)) 
— ^ — 



-i 



(43) 



Using (l4"3l . we can derive the slope factor A, 



(44) 

which can be used to estimate the variance using equation ( [211 . 

We also evaluate the response function numerically. We assume that the response function 
can be approximated by the firing probability of a neuron whose membrane potential obeys 

where X G {0, ...,iV}. This presumes that the effect of discrete synaptic events can be 
approximated by a constant input equivalent to the mean input plus noise. We estimated 
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the firing probability in a temporal epoch from the average firing rate of the neuron obeying 
(|45p for each X. We used the firing probabilities directly to form the MTM. We then 
computed the equilibrium PDF, the mean activity (crossing point), the slope factor, and 
the variance. We choose the bin size of the Markov model to be the membrane time constant. 
All differential equations were computed using the Euler method with a time step of 0.01 
ms. 

We compared the numerically simulated mean and variance of the network activity at 
equilibrium for 100 neurons with our theoretical predictions for differing synaptic strength 
J. The results are displayed in figure [HJA A close match is observed for both quantities. 
The mean and variance increase with J. Figure [HJA. also shows the mean- field (uncorrelated) 
variance. We see that the network variance always exceeds the mean-field value especially 
for midrange values of J. The PDFs from the numerical simulation of 100 neurons and using 
the Markov model are shown in[HjB. The simulated PDF was generated by taking a histogram 
of the network activity over 10 5 ms. For the Markov model, the PDF is the eigenvector with 
eigenvalue one and approximated by taking a single row of the one hundredth power of the 
MTM. 

The model predicts an exponentially decaying autocorrelation function with time constant 
l/ln(A). It is displayed in figure [HJC for a network of N = 1000. We used the same 
parameters as in figure [HA with a connection weight where the dynamics deviates from 
mean-field theory (J = 5). In this parameter range, the recurrent excitation is strong and 
the single neuron firing rate is very high, approximately 350 Hz (the refractory time of 1 
ms imposes a maximum frequency is 1000 Hz). In this regime, the refractory time acts 
like inhibition, so the probability to fire actually decreases if the activity in the previous 
epoch increases. Thus, the autocorrelation function exhibits anti-correlations as seen in the 
figure [HJC. We can estimate A by fitting the autocorrelation to A T . Using the approximation 
A = Cov(l)/Var, which gives |A| = 0.40. The theoretical value of |A| using f j44l) gives 0.37 
for the same parameters (using the simulated mean number of firing neurons). Figure 03 
shows a plot of the variance versus N for the numerical simulation. The simulated variance 
is well matched by the estimated variance from (12T1) with q measured from the simulation 
(estimated by the mean firing rate of a single neuron) and A = 0.4. 
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Conductance-based synapse 



We now consider the integrate-and-fire neuron model with conductance-based synaptic 
connections. The membrane potential V obeys 

dV 

r— = I(t) -{V- V r ) - s(t) (V - V.) + Z(t) (46) 

where r is the membrane time constant, I(t) is an input current, Z(t) is zero-mean white 
noise with variance a, V r is the rest potential, and V s is the reversal potential of an excitatory 
synapse. The synaptic gating variable s(t) obeys 

r. d i=±p(t-t')s(t) (47) 

where t s is the synaptic time constant, J is the maximal synaptic conductance, N is the 
number of neurons in the network, and v are the firing times of all neurons in the network. 
The threshold is Vg and the reset potential is Vr. As with the previous model, a refractory 
period r was introduced. We used r = 1 ms, t s = 1 ms, V r = —65 mV, V s = mV, Vg = —60 
m V, Vr = —80 mV, t s = 1 ms and r = 1 ms. 

We computed the response function by measuring numerically the firing probability of a 
neuron that obeys 

= I ~ (V - V r ) ~^(V- V rvs ) + Z(t) (48) 

for all X G {0,...,iV} using the same method as in the current-based model. From the 
response function, we obtained the MTM and the invariant measure. Figure compares 
the mean and variance of a numerical simulation with the Markov prediction for varying J. 
We see that there is a good match. For J = 0.05, we compare the numerically computed PDF 
with the invariant measure. As shown in figure [9j3, the prediction is extremely accurate. 
We estimated the slope factor as in the previous section using the autocorrelation function 
for iV = 1500 and found A = 0.31 and computed the estimated variance for various N. The 
comparison is shown on figure [9p. There is again a close match. 

DISCUSSION 

Our model relies on the assumption that the neuronal dynamics can be represented by a 
Markov process. We partition time into discrete epochs and the probability that a neuron will 
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fire in one epoch only depends on the activity of the previous epoch. For this approximation 
to be valid, at least two conditions must be met. The first is that the epoch time needs to be 
long enough so that the influence of activity in the distant past is negligible. The second is 
that the epoch time needs to be short enough so that a neuron fires at most once within an 
epoch. The influence time of a neuron is given approximately by the larger of the membrane 
time constant (r m ) and synaptic time constant (r s ). Presumably, the neuron does not have 
a memory of events much beyond this time scale. This gives a lower bound on the epoch 
time At. Hence, At > max[r m ,r s ]. The second condition is equivalent to / At < 1, where / 
is the max firing rate of a neuron in the network. Thus, the Markov model is applicable for 

max[r m , r 8 ] < f' 1 (49) 

This condition is not too difficult to satisfy. For example, a cortical neuron receiving AMPA 
inputs with a membrane time constant less than 20 ms can satisfy the Markov condition if 
its maximum firing rate is below 50 Hz. Since typical cortical neurons have a membrane 
time constant between 1 and 20 ms and many recordings in the cortex find rates below 50 
Hz 40J, our formalism could be applicable to a broad class of cortical circuits. 

The equilibrium state of our Markov model exists and is exactly solvable if the response 
function is never zero or one. In other words, a neuron in the network always has a nonzero 
probability to fire but never has full certainty it will fire. Hence, the neuronal dynamics 
are always fully stochastic. Thus the equilibrium state of a fully stochastic network could 
be another definition of the asynchronous state. However, even though the network is 
completely stochastic, the activity is not uncorrelated. These correlations are manifested 
in the entire network activity but not within the firing statistics of the individual neurons, 
which obey a simple random binomial or Poisson distribution. The autocorrelation function 
of the individual neuron also decays much more quickly than that of the network. 

Many previous approaches to studying the asynchronous state assumed that neuronal 
firing was statistically independent so a mean-field description was valid. With our simple 
model, we can compute the equilibrium PDF directly and explicitly determine the parameter 
regimes where mean-field theory breaks down. We also show that the "order parameter" that 
determines nearness to criticality is the slope of the response function around the equilibrium 
mean activity. As expected, we find that mean-field theory is less applicable for small or 
strongly recurrent neural circuits. In our model, the mean activity of our network can be 
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obtained using mean-field theory except at the critical point. We compute the variance 
of the network activity directly and see precisely how the network dynamics deviate from 
mean-field theory as we approach the critical point. We note that while our closed-form 
estimates for the mean and variance may break down very near criticality, our model does 
not. The invariant measure of the MTM still gives the equilibrium PDF and in principle 
can be computed to arbitrary precision. Hence, our model can serve as a simple example of 
when mean-field theory is applicable. We note that even in the limit of N going to infinity, 
the variance of the network deviates from a network of independent neurons by a factor of 
(1 — A) -1 . The deviations from mean-field theory persist even in an infinite sized network. 

Our model requires statistical homogeneity. While purely homogeneous networks are 
unrealistic for biological networks, statistical homogeneity may be more plausible. For some 
classes of inputs, the probability of firing for a given set of neurons could be uniform over 
some limited temporal duration. We found that our analytical estimates can be modified 
to account for disordered randomness in the connections. A model that fully incorporated 
all heterogeneous effects would require taking into account different populations. However, 
with each additional population, the dimension of the MTM increases by a power of N. 
Thus for a network of two populations, say inhibitory and excitatory neurons, the resulting 
problem will involve a MTM with iV 2 x iV 2 elements. 

Future work could examine the behavior at the critical point more carefully. Our variance 
estimate predicted that the variance at criticality would diverge as the system size squared 
but simulations found an exponent of 3/2. We also showed that the correlation time of 
the network activity diverged at the critical point. We expect correlations to obey a power 
law at criticality. Perhaps, a renormalization group approach may be adapted to study the 
scaling behavior near mticality. Recent experiments have found that cortical slices exhibit 
critical behavior 
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42]. Our Markov model may be an ideal system to explore critical 



behavior. 

Finally, we note that even at the critical point, where the network activity is highly 
correlated, the single neuron dynamics can exhibit uncorrelated Poisson statistics. This 
shows that collective network dynamics may not be deducible from single neuron dynamics. 
Using our model as a basis, it may be possible to probe characteristics about local network 
size, connectivity and nearness to criticality by combining recordings of single neurons with 
measurements of the local field potential. 
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MEAN AND VARIANCE DERIVATIONS 



The mean activity is given by 

N 

(X) t = J2 kP k (t) 

k=0 

N N 



N N 

N-kk 



k=0 8=0 

Rearranging yields 

(*)t = E^-i)E kC N(l - Pi) N ~ k p 

i=0 k=0 

Since Y.k=okC^(l — Pi) N ~ h Pi = Npi is the mean of a binomial distribution, we obtain 

N 

(X) t = Y,Pi(t-l)Npi 
= N(p) t ^ 

which is (JHD- 

The variance is given by 

N 



Var t (X) = Y.k 2 P k (t)-N 2 (p)U 

k=0 
N N 

= E fc2 E c h N (i - Pi ) N - k P k m - 1) - n^ p )U 



k=0 i=0 
N N 

= E m - 1) E k 2 c k N {i - Pl f- k P k - n\p)U 

i=0 fe=0 
N N 

= E m - i)N Pi {i -p*) + E Pit - l ) N Vi - n 2 (p)U 

i=0 i=0 

= N{p{l-p)) t _ l + N 2 Var t _ 1 {p) 

where we have used the variance of a binomial distribution Np(l — p). For the linear case, 
writing p(X) = a + bX we have 

Var^X) = N(p(l-p))^ + N 2 (p 2 )^N 2 (p)l 

= N(a + bX - a 2 - 2abX - b 2 X 2 )^ + N 2 (a 2 + 2abX + b 2 X 2 )^ -N 2 (a + 6X) 2 
= Na + bN(X) - No 2 - 2Nab(X) - Nb 2 (X 2 ) + N 2 a 2 + 2N 2 ab(X) + N 2 b 2 (X 2 ) 

-N 2 (a 2 + 2ab(X) +b 2 (X) 2 ) 
= Na- Na 2 + N 2 a 2 - N 2 a 2 + (X)(bN - 2Nab + 2N 2 ab - 2N 2 ab) 
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+ (X 2 )(-Nb 2 + N 2 b 2 ) - b 2 N 2 (X) 2 

Na - Na 2 + bN(l - 2a) (X) + (N 2 b 2 - Nb 2 )Var(X) - Nb 2 (X) 2 
Na - Na 2 + bN(l - 2a) (X) - Nb 2 (X) 2 
1 - N 2 b 2 + Nb 2 



Setting a = p and b = ^ A/ po ' ) gives 



^° r " (X) = 9 W «"-(«-«,)' + («-«,)>/" (50) 



If we set A = (g — po)/g we get equation (12T]) . 



AUTOCOVARIANCE FUNCTION 

We prove the form of the autocovariance function Cov(t) = X T Var fl (X) for the linear 
response function using induction. We first show that Cov(l) = XVar^(X) and then Cov(r+ 
1) = \Cov{t). 

The autocovariance function when r = 1 is given by 

Cov(l) = (X t X t+1 ) - (X); 



N N 

2 
I" 



J2 E k ^P( x t+i = i\X t = k)P(X t = k) - (X) 

k=0 i=0 



= Nj2k Pk P(X t = k)-(X)l 

k=0 

= N(pX), - (X)l 
For a linear response function p, we obtain 

Cov(l) = N(pX), - (X)l 

= N( P0 X + q —^X 2 ), - N 2 m 2 
JSq 

at-2 . /Q ~~ P0 v 2\ 7\r2 2 

= N p q + ( — - — X ) lt -N q 



N 2 q (po - ?) + ^— ^ (Var M (X) + iVV 



Var M (X) 

Hence for r = 1, the autocovariance function is equal to the slope factor X — (q — po)/q 
times the variance. Assume for r, that Cov(t) = X T Var ll {X), then 

Cov(r + l) = (X t X t+T+1 ) - (X)l 
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N N 

= E E k J p ( x t = k)P(X t+T+1 = j\X t — k) — (x)l 

k=0 j=0 
N N N 

= E E k J P ( X t = k ) E W+r+i = j\X t+T = r)P(X t+T = r\X t = k)- (X)l 

k=Oj=0 r=0 
N N N 

= E E E = - p(r)) JV -^(r)^P(X t+T = r\X t = k) - {X% 

k=0j=0r=0 
N N 

= E E Nkp(r)P(X t = k)P(X t+T — r\X t = k) - (X)J 

fc=Or=0 

using the mean of the binomial distribution. We can now insert the linear response function 
to obtain 

^ N a - d 

Cov(t + 1) = E E Nk(po + \- ^r)P(X t = k)P(X t+T = r\X t = k) - (X) J 

N N 

= £ £ XA;p P(X t = A;)P(X t+T = r|X t = *) + *— ?S (CW(r) + (X) J) - (X); 

fc=0 r=0 9 

because 

iV 

\2 



E E fcrP (^ = k)P(X t+T = r\X t = k) = Cov(t + r) + (X) 

k=0 r=0 



Then since 



N N N 

E E Nk Po P(X t = k)P{X t+T = r\X t = k) = Y, Nk Po P(X t = k) 

k=0 r=0 k=0 

because 

N 

Y / P(Xt+r = r\X t = k) = l 

r=0 

for all k and 

N 

E Nkp P(X t = k) = (Np X)^ = N 2 p q 

k=0 



we finally obtain 



Cov(t + 1) = N 2 p q + ?—^(X)l-(X)l + \Cov(T) 

= N i poq + 1^N 2 q 2 - N 2 q 2 + XCov(r) 
Q 

= N 2 (p q + (q- p )q - q 2 ) + XCov(r) 



= XCov(t) 
proving equation ff24l) by induction. 
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Figures 




FIG. 1: Example for response function p(n) = ^- fe^i-j n e ? dx with N = 100, 9 = 1, I = 

a 

0.1, J = 1.5/N and a = 1.0. A) Np(n) (solid line) and the diagonal (dashed line). B) The 
associated Markov Transition Matrix (in gray levels). C) The evolution of the mean network 
activity (solid line) and the variance (dashed line). The steady-state is quickly reached and the 
activity corresponds to the crossing point between Np(n) and the diagonal. D) The invariant 
measure (i.e the PDF of the network activity). 
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FIG. 2: Example where the response function crosses the diagonal three times. A) Np(n) (solid 
line) and diagonal (dashed line). B) The associated MTM. C) The invariant measure. D) MTM 1q5 . 
All rows are equal to the invariant measure. Parameters are N = 100, I = 0.1, J = 1.8/N and 
a = 0.6. 
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FIG. 3: A) The two dimensional bifurcation diagram (A = 1) as a function of and — . Inset: 
the traversal lines of the various one dimensional bifurcation diagrams shown in B) vertical, C) 
horizontal and D) oblique. 
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FIG. 4: A) Numerical and theoretical PDF of the network activity at equilibrium for N = 100, 
8 = 1, I = 0.1, a = 0.8 and J = 1.8. The theoretical PDF was obtained by taking one row of 
MTM 100 . B) Numerical and theoretical autocorrelation function for same parameters as A). Solid 
line is A* with the predicted A from equation (|30H . C) Dependence of the mean (solid line) and 
variance (dashed line) with J. The mean-field variance is Nq(l — q) (dotted line). D) Evolution 
of the variance when A = 1 versus N for 1 = 1 — V27T, J = V2tt, o = 1, 9 = 1. The solid line is 
N 3 / 2 q(l-q). 
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FIG. 5: Autocorrelation of the network and one neuron at a critical point (A = 1). The solid line 
is for a network of 10,000 neurons and the dashed line is for a network of 1000 neurons. Crosses 
and dots are the autocorrelation of one neuron for N = 1000 and N = 10, 000 respectively. The 
autocorrelation for one neuron decays almost immediately. Parameters are a = 1, 9 = 1, J = \p2/n 
and 1=1 — ^/2j^. 



29 




1000 2000 3000 4000 5000 10 20 30 40 50 

Time Number of Firing Neurons 




200 400 600 800 1000 20 40 60 SO 100 120 140 

Number of firing neurons Number of neurons 



FIG. 6: A) The firing activity over time of a bistable network of 50 neurons. The switching 
between the two states (high and low activity) is spontaneous. B) Theoretical PDF for N = 50. 
C) Theoretical PDF for N = 1000. D) Switching rate depends exponentially on ./V (solid curve). 
Dashed line is a linear fit on semi-log scale. Parameters are: J = 1.8, a = 0.6, I = 0.1, 9 = 1. 
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FIG. 7: Mean (crosses) and variance (diamonds) from a numerical simulation of the fast leak model 
for N = 100 neurons as a function of the variance of the random disordered connection weights. 
Theoretical mean (solid line) and variance (dashed line) match well with numerical simulation 
values. Parameters are J = 0.5/iV, / = 0.1, a = 1.0, 9 = 1.0. 
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FIG. 8: Current Synapse model - A) Mean (crosses) and variance (diamonds) versus J from the 
numerical simulation compare well with the predictions (solid and dashed lines). The mean- field 
variance is also plotted (dotted line). Parameters are N = 100, I = 4, a = 2.0. B) The numerical 
PDF(histogram) and prediction (solid line) for N = 100, J = 5.0, 1 = 4, and a = 2.0. C) The 
autocorrelation (solid line) for a network of N = 1000. The dashed line is the estimated exponential 
decay with A = Cov(l)/Var (see text). Parameters are J = 5, I = 4, a = 2.0 and A = 0.40. D) 
The variance (circles) versus N. The solid line is from equation (|2ip with A computed using the 
autocorrelation decay for iV = 1000 (A = 0.4) displayed in C). Parameters are J = 5, I = 4, 
a = 2.0. 
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FIG. 9: Conductance Synapse model - A) Mean (crosses) and variance (diamonds) versus J from 
the numerical simulation compare well with the predictions (solid and dashed lines). Parameters 
are N = 100, I = 4, a = 2.0. B) The numerical PDF(histogram) and prediction (solid line) for 
N = 100, I = 4, a = 2.0, J = 0.05. C) The autocorrelation (solid line) for a network of N = 1500. 
The dashed line is the estimated exponential decay with A = Cov(l)/Var (see text). Parameters 
are J = 0.05, I = 4, a = 2.0 and A = 0.31. D) The variance (circles) versus N. The solid line 
is from equation (|21|) with A computed using the autocorrelation decay for ./V = 1500 (A = 0.31) 
displayed in C). Parameters are J = 0.05, I = 4, a = 2.0. 
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